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Abstract 

In this paper we study the on and off-equilibrium properties of the four dimen- 
sional Gaussian spin glass. 

In the static case we determine the critical temperature and the critical exponents 
with more precision that in previous simulations. 

In the off-equilibrium case we settle the general form of the autocorrelation func- 
tion, and show, for the first time, that is possible to obtain, dynamically, a value for 
the order parameter. 



1 Introduction 



Nowadays, the problem of the full characterization of the phase transition in finite dimen- 
sional spin glasses is still opened both from the static and dynamical approaches. 

Our discussion is focused on the four dimensional case (the same applies in the more 
physical case, the three dimensional system). 

The equilibrium (static) simulations show a very neat intersection of the Binder cu- 
mulant curves, that is a signal of a phase transition at finite temperature with an order 
parameter (the Edward- Anderson order parameter, that we will denote hereafter as qEA)- 
We can identify this order parameter with the position of a Dirac delta in the probability 
distribution of the overlap, P(q). Up to now, both the SRSB theory (Spontaneous Replica 
Symmetry Breaking) |T], |2|] and the droplet theory || are compatible with this result. Dif- 
ferences concern the shape of the rest of the P(q). In the droplet theory P(q) is the sum 
of two Dirac deltas, one in q EA and another in the opposite overlap, and has a Binder cu- 
mulant equal to 1. The SRSB theory maintains this structure too, but adds a continuous 
non zero part in the interval (— qEA, qEA)- This is a non trivial distribution that has a 
Binder cumulant different from 1, except at T = where the SRSB theory predicts two 
pure states like the droplet theory. 

The main problem is the impossibility of a direct measure of the order parameter, 
qEA- The scaling of the peak of P(q) seems compatible both with a Kosterlitz-Thouless 
(KT) transition, i.e. q pea k ~ 1/L a , and with a scaling like q pea k = qEA + a/L p . Obviously 
the KT scenario goes against the intersection of the Binder cumulant curves. A possible 
explanation of this phenomena could be that the term a/L p is bigger than qEA for the range 
of lattice sizes that has been simulated and hence the latter is unobservable. Simulation 
of bigger lattices should be done in order to get qEA S> a/L p . 

The dynamical approach M has the same problematic of the previously discussed static 
case. The main physical quantity in this approach is the spin-spin autocorrelation, defined 
as: 

1 N 

C{t,t w ) = -Y,(vi(t w )vi{t w + t)) , (1) 

iV i=l 

Usually in the literature || one finds the empirical formula 

C(t,t w ) =t- x f(t/t w ) , (2) 

for instance in the Mean Field case and in the three dimensional case 0Q. 

The static limit (on-equilibrium situation) is achieved sending t w to infinite first, and 
then simulating larger t. The formula (0) for the spin-spin autocorrelation function goes to 
zero in this limit. However, in the case of a non zero order parameter, this autocorrelation 
function must go to qEA- It is clear that in the regime of t w t ^> 1 should be found a 
formula like 

C(t, t w ) = (q EA /f(0) + ar x )f(t/t w ) , (3) 

but a very long numerical simulation is needed in order to observe both the terms q E A and 
at~ x . In the present work we show for the first time numerical evidence of this kind of 
behaviour. 

1 Also in an on-equilibrium numerical simulation in the three dimensional case C(t) ~ t~ x . 
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Up to now, the only numerical studies of off-equilibrium dynamics in finite dimensional 
spin-glass are those of H. Rieger in the three dimensional case. 

The four dimensional case seems easier to simulate, since it is far away from the lower 
critical dimension of the spin glasses (di < 3 || [TO]) and the static is thus very clear. In this 
paper we will study mainly the off-equilibrium dynamics of this model in order to compare 
with the three dimensional results by Rieger and to examine the possibility to extract a 
finite value for the order parameter. In addition simulations have been performed in the 
static (on-equilibrium) case in order to characterize with higher precision the location of 
the critical temperature and the critical exponents. 

Both for the off and on-equilibrium cases we review the numerical results from the point 
of view of the previous discussion and we try to link both approaches in order to obtain 
a conclusion regarding the existence of a finite temperature phase transition, with a non 
zero order parameter. 



2 Model, simulation and static observables 

We have studied the 4-d Ising spin glass with nearest neighbor interactions and zero external 
magnetic field, whose Hamiltonian is 

H = - ■ l 'j a ' a .i > ( 4 ) 

<i,j> 

where < i,j > denotes nearest neighbor pairs and the couplings are extracted from a 
Gaussian distribution with zero mean and unit variance. 

The static and dynamical behaviour of the model have been investigated by several 
different simulations during which many different observables have been measured . This 
section describes the way we performed measurements of the static exponents and the 



critical temperature with a precision higher than that available in the literature pl| , |12[ . 

The equilibrium simulations have been performed on small lattices (L = 3, 4, 5, 6, 7, 8) 
to ensure the system reached equilibrium. Most of the work has been made in a range of 
temperature around the critical one, T c . The average over the disorder, has been done on 
2048 samples for all the lattice sizes. For each realization of the quenched disorder we have 
simulated two replicas with spin (jj and Tj. This enabled us to measure the k-th cumulant 
of the distribution of the overlaps, = J q h P(q)dq, simply by averaging the quantity 
(N^ 1 J2i o r iTi) h over a large number of independent configurations. 

All the calculation have been carried on a tower of the parallel supercomputer APE100 
fill , with a real performance of about 5 Gigaflops. 

A detailed study has been devoted to the calculation of the number of sweeps needed to 
reach the equilibrium and to the estimate of the autocorrelation time at the equilibrium. 
This study suggests a thermalization time of about 10 5 sweeps, being sure that using this 
value even the biggest system at the lowest temperature will be thermalized. To verify 
the correctness of this value we studied the evolution of the biggest system (L = 8) at the 
lowest temperature (T = 1.7): we choose three replicas of the system such that having, at 
the starting time, two overlaps set to zero and the third one equal to one; we have followed 
the evolution of these overlaps averaging over a large number of disorder configurations 
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Figure 1: xsg{L,T) vs. L\ the errors are of the order of the symbol. 



and we have estimated the thermalization time as the time needed in order that the three 
overlaps converge to a single value. 

Once the equilibrium has been reached, we measured how much time was needed to 
decorrelate the observables. Particularly, we have seen that the overlap between two repli- 
cas has a time correlation function that decreases exponentially, C(t) ~ exp(— t/r). This 
defines a characteristic time whose typical value at T = 1.7 are: for L = 4 r ~ 200, for 
L = 6 t ~ 1000 and for L = 8 r ~ 3000. In the final simulation, after thermalization, 
we measured every r sweeps the overlap between the two replicas for a time longer than 
the equilibration one. 

Defining the spin glass susceptibility as 



1 



XsG (L,T) = -J2(^) 2 = Nq 

i,3 



(2) 
L 



(5) 



where N = L 4 , ((••)) is the thermodynamical average and (••) the mean over the disorder; 
and the Binder parameter as 



1 / o (4) 
g(L,T) = - 3-4^ 

2 I (g « )2 



(6) 



the results of our simulations are plotted in fig.|l] for the spin glass susceptibility and fig.0 
for the Binder cumulant. 

The errors on the plotted data are derived from a jackknife analysis, which also confirms 
that the overlaps measured every r sweeps are decorrelated. Using finite size scaling we 
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Figure 2: g(L,T) against T; the errors are of the order of the symbol. 



see that xsg(L,T) and g(L,T) scale as (in the scaling region) 

X sg(L,T) = L 2 ^xsg(L 1/v (T - T c )) 
g(L,T)=g(LV»(T-T c )) . 



(7) 



Note that at the critical temperature the Binder parameter does not depend on the size of 
the system, so T c is the temperature where the curves of fig.0 cross themselves. 

In the neighborhood of T c we can approximate the function g with a linear one and 
obtain the following critical temperature and v exponent 



T 

-L C 

V 



1.80 ±0.01 
0.9 ±0.1 . 



(9) 
(10) 



The value of v is also confirmed by the results of the analysis done, following fl4}| , on the 
quantity 

dg 



dT 



obtaining 



To: g(T )=g 

1.06 ±0.06 



(11) 
(12) 



The prediction about the infinite volume limit of the Binder cumulant is different in the 

™g(T)<i). 

3 to L = 



droplet theory (g(L, T < T c ) Z ^? > 1) and in the SRSB picture (g(L,T <T C ) - 



Unfortunately with our data (L = 3 to L = 8) is impossible to extrapolate the infinite 
value with a precision that could discriminate, without doubt, between the two predictions. 
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The estimation of the anomalous dimension 77 can be performed by fitting the xsg data, 
at T = T c , with a power law 



Xsg(L,T = T c ) ocL 



2-7) 



(13) 



obtaining rj = —0.35 ± 0.05 (the error is due mostly to the indetermination on the crit- 
ical temperature and to the rapid variation in the region around T c of the exponent in 
eq . (|i~3l) ) . These results are in agreement with those found by Bhatt and Young in [11 
using a maximum size of 6 4 and 200-800 samples: T c = 1.75 ± 0.05, v = 0.8 ± 0.15 and 
rj = -0.3 ±0.15. 

Using the scaling law 7 = v{2 — rj) and the exponents values just calculated, we have 
7 = 2.1 ±0.2, which is in good agreement with the value obtained by the high temperature 
expansions 7 = 2.0 ± 0.4, |T5|. 



Another series of computer runs, performed using the annealing procedure let us 
measure the non-connected susceptibility for a wide range of temperature in the spin glass 
phase (T < T c ). We clearly see that the data diverge with increasing system sizes, even 
though, because of the small lattices, many different fits are possible, e.g. 



X5G (L,T) = A(T)L 4 1 +B(T)L 



-A(T) 



(14) 



or 



Xsg{L,T) oc L 



2-n(T) 



(15) 



Further evidence for the value of T c can be obtained, as the highest temperature where the 
power law fit is yet acceptable (by a \ 2 test). 



3 Off-equilibrium dynamics 

The second part of our study was devotes to the simulation of systems of greater dimensions 
pranging from 8 4 to 32 x 16 3 ). At the be ginning of every simulation the system is frozen 
from an infinite temperature to one in the critical region (T < T c ), and measures of the 
autocorrelation functions immediately start with the system still out of equilibrium. Due 
to the huges thermalization times of the simulated systems, the off-equilibrium dynamics 
is the most realistic situation and also the most interesting. In fact, due to the enormous 
number of metastable states, the dynamics is very slow and besides it is reminiscent of the 
time passed in the spin glass phase, that we will call t w . These effects can be clearly seen 
by the study of the autocorrelation functions defined as 

1 N 

C(t,t w ) = —J2(ai(t w )ai(t w + t)), (16) 

where (••) is the mean over the disorder and ((••)) stands not for an average over the 
equilibrium thermodynamic state, since we are not at the equilibrium, but for an average 
over the thermal histories. Nevertheless we found that, for the system sizes we considered, 
disorder fluctuations are always stronger, so generally we omit the angular brackets. 
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Figure 3: C(t, t w ) vs. t at T = 0.2 with t w = 2 7 , 2 10 , 2 13 , 2 16 , 2 19 (bottom to top). 

Our simulations cover the cold phase (from T = T c = 1.8 down to T = 0.2) through 
the set of waiting times t w = 2 k with k = 7, 8, . . . , 21 and averaging over 3072 disorder 
realizations systems of volumes from 8 4 to 12 4 . 

In the four-dimensional Ising spin glass the presence of a critical temperature and the 
subsequent spin glass phase has been widely accepted, so the principal question that re- 
main to answer is which kind of phase space arises for T < T c . In the literature there are 
principally two theories that try to describe the spin glass systems in their low tempera- 
ture phase: one is based on a mean field like approximation which predicts a spontaneous 
replica symmetry breaking (SRSB picture); the other one, starting from a Migdal-Kadanoff 
renormalization group technique, concludes that the system remains trivial, with only one 
pure state (droplet model). The predictions of the two theories regarding the autocorre- 
lation function are different: the SRSB picture predicts that in the limit of t w — > oo the 
autocorrelation must be a power law that converge to the Edward-Anderson parameter 
(Qea) 

C(t,t w ) = (q EA + ar x )^jl^ , (17) 
while in the droplet model the relaxation is slower 

„/l0g(t/T)\ 



C(t,t w ) = (\ogt)~ e ^ c 



log(W r ) 



The data we collected (see fig|| and fig.|]) seem to agree with the scaling law used in 
and in 

C(t,t w )=t- x '^C(t/t w ) , (19) 
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Figure 4: C(t,t v 



vs. 



t at T = 0.45 with t w = 2 7 , 2 10 , 2 13 , 2 16 , 2 19 (bottom to top). 



with a scaling function 



C(z) 



constant 

z x'(T)-\(T) 



for z 
for z 





oo 



(20) 



The values of the exponent x'(T) is plotted together with the values of x(T) in fig.0, while 
A(T) is plotted in fig.[7| . 

To evaluate the goodness of the two proposed scaling formulas eq.(|19"D and eq.flT8|) we 
plotted in fig.[| the T = 0.45 data rescaled with the former law, noting that they collapse 
very well on a single curve. On the contrary, using the droplet model scaling law, it was 
impossible for us to find a value for the parameters 6/ip and r such to force the data over 
a single curve; in fig.|^ we show the best scalings we could obtain in order to make the data 
collapse in the t < t w or in the t > t w region. 

Nevertheless the very good rescaling of the data in fig.|5] we also performed a deeper 
analysis in order to find the value of the Edward- Anderson parameter, qEA-, which is as- 
sumed to be zero in eq. (|I9"D . The value of qEA can be found performing the t — > oo limit 
after the t w — > oo limit; for this purpose we have done very long simulations (more than 
4 millions Monte Carlo sweeps). We note that the scaling laws obeyed by the data in the 
two regions t ^> t w and t <C t w are essentially different. In the former the data can be 
fitted by a power law of the ratio t/t w , while in the latter we obtain a law equal to that of 
eq.(|T7D times a function of t/t w which is almost a constant. 

Such a behaviour for C(t, t w ) can be justified assuming that the system evolves as long 
as t <^ t w with a quasi equilibrium dynamics which converges to qEA while for t 3> t w it 
decorrelates faster and toward zero (C ~ t~ x with A(T) ^> x(T) VT), but always with a 
critical slowing-down. 
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Figure 5: T = 0.45 aging autocorrelation function rescaled following eq. (|19|) , with 
x = 0.0054. We plot t x C(t,t w ) against t/t w . 
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Figure 6: Tentatives of rescaling following the droplet model law eq. (|i8|) : in the left plot 
e/ij) = 0.0054 and log(r) = -1; in the right plot 6/tp = 0.043 and log(V) = -1000. 
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Figure 7: X tw =o(T) (top) and Xt w ^o(T) (bottom) against T; the lower line represents the 
power fit A(T) = 0.21(1) T 0,92 ^; the upper line both the linear and the power fit: respec- 
tively A(T) = 0.000(3) + 0.30(1) T and A(T) = 0.303(8) T im ^. 

The values for A(T) have been obtained fitting the C(t,t w = 0) data with a power law 
and for t w ^ with 

/ t \ ~ A ( T ) 

C(t,t w )<x - , (21) 

in the range t/t w > 15. In fig.[7| we plot the results either for t w = and for t w ^ 0. We 
note that both fits are compatible with the linear dependence in the temperature predicted 
from the experimental measurements ||17|| . 

In the region t w /t > 32 we have performed the analysis assuming that the correlation 
function could be factorized as 

C(t,t w ) = (q EA + at~ x )C(t/t w ) , (22) 

where we have approximated C(z) = l — ciz C2 for z — > 0. First of all, the rescaling function 
C(t/t w ) has been calculated fitting the correlation function at a fixed value of t. Later, 
once divided the data by this function, we verified that the curves for differents ratios t/t w 
collapse over a single curve and we interpolated the data via a power law plus constant, 
following eq. (p~7|) . In fig|| we plot in a log-log scale typical C(t,t w )/C(t/t w ) data with the 
best fit; we note that up to now in the literature this data have been fitted via a simple 
power law, while it is evident that the points in fig.| are not on a straight line. 

From these fits the values of q E A and i as a function of the temperature (see fig.[5] and 
fig.|K]) can be obtained. As a guide to the eye, we plot in fig.pl] the simpler function that 
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data at T = 0.9 versus t; the upper line is the best power plus constant 



fit: 0.60(4) + 0.32(4) £-°- 08 W, while the lower line is the best power fit. 



behaves like \T — T c f near the critical temperature and tends to 1 for T = 



Qea(T) 



(23) 



where T c = 1.8 and (3 = |(c2 — 2 + 77) = 0.74 (using the values found in the previous 
section). From fig|| we note that only the quantity x(T), and not x'(T), is such that 
x(T)/T is roughly independent from the temperature, so that only in this parametrization 
the t w = 00 autocorrelation function (R(t;T) = C(t,t w = 00) at temperature T) can be 
written as 

R(t;T) - R(oo;T) = b(T)exp(-BT\og(t)) . (24) 

The relevance of the variable Tlog(t) has been observed in experiments on magnetic re- 
manence in a wide region |17| . 

We call off-equilibrium correlation length, £(£), the typical distance over which the 
system is thermalized after a time t. For this domain growth the SRSB picture predicts a 
power law || 

m oc t 1/z(T) , (25) 

while in the droplet model, where the energy barriers scale proportionally to Lr, the law 
is 

^(t)oc(Tlogt) 1 ^ . (26) 
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Figure 9: x(T) (top) and x'(T) (bottom) against T (the value at the greater temperature 
is equal: x(T c ) = x'(T c )); the line is best power fit x(T) = 0.083(3) T im<y7 \ 




Figure 10: Edward- Anderson order parameter against T; the line is only a guide to the eye 
as explained in the text. 
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At least at the critical temperature there is a scaling relationship between the dynamical 
exponent z(T c ) and the one which describes the dynamics in the quasi equilibrium regime, 
x(T c ) = x'(T c ) (because q E A(T = T c ) = 0) 

This equation is satisfied by all the exponents we have estimated at the critical temperature: 
x = 0.15, i] = —0.35 and z = 5.3. 

To find the behaviour of the off-equilibrium correlation length we have measured, like 
in ||, the equal time spatial correlation functions 

1 N 



G(r,t) = -J2W l (t)a l+r (t)) 2 , (28) 

ly i=l 

where the averages are the same as in eq.(|l^) and t is the time since the cooling. This 
study has been performed on systems of volume 32 x 16 3 . 

From scaling concepts we know that, at large values of r, G(r,t) must behave like 

G(r,t)ocr-^- 2+ ")/^^ , (29) 

and, supposing f(y) = Aexp(—By D ), we have fitted our data with the function 



G(r, t) = Ar- {2+Tl) 



exp 



-B 



r 



(30) 



In the t — > oo limit the exponential term tends to 1 and we obtain a spatial correlation 
function that decrease with a power law: in fig.|ll] we plot such function at the critical 
temperature (T c = 1.8). Note that from the slope of the curve we obtain an estimation of 
the r) exponent compatible with that of section ^|. 

At the lower temperatures the value of r\ strongly depends on the r range of interpola- 
tion, because the fitting function diverges at r = 0. On the contrary, trying to fit the data 
in different ranges of r, we find that the dynamical exponent z(T) is a robust parameter 
which remains unchanged for every r range (plotted in fig.0) 

Fitting the plotted data with a power law we obtain, up to the critical temperature 

z(T) = A T~ a , (31) 

with A = 9.7 ± 0.5 and a = 1.0 ± 0.1. 

A preliminary analysis of a new set of data at T = 0.9 = T c /2 suggests a value of 
rj ~ — 1. The fact that the value of z at this temperature is higher than the corresponding 
value at T c makes the evaluation of the exponent rj prone to systematic error. Nevertheless 



this rough estimate of rj is compatible with the prediction of the reference |L8 . 
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Figure 11: G(r,t = oo) vs. r at T = T c = 1.8; the line is the power fit 0.24(1) j- -1 - 63 ®. 
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Figure 12: z(T) vs. T; the line is the power fit, eq.( 
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4 Conclusions 



In this paper we have studied the on and off-equilibrium properties of the four dimensional 
Gaussian spin glass. 

In the static case the hypothesis of a transition a la Kosterlitz-Thouless have been 
rejected owing to the study of the on-equilibrium properties of the model: the existence of 
a second order phase transition is well testified by the clean cut of the Binder cumulant 
curves. We have determined with more precision that in previous simulations both the 
critical temperature as well as the critical exponents. 

In the off-equilibrium case we have settled, for the first time, a form of the autocorrela- 
tion function compatible, in the large times limit {i.e. on equilibrium), with the existence 
of an order parameter different from zero. We have been able to determine, in a dynamical 
way, the value of q^A as a function of the temperature (this again confirm the absence of a 
Kosterlitz-Thouless transition). Also we have established the temperature dependence of 
the exponents that appear in these off-equilibrium dynamics, linear in all the cases. 

The dynamics of the model seem to be described better by the SRSB theory than by 
the droplet theory: in effect the autocorrelation functions, properly rescaled, follow very 
well the power laws predicted by the former, while are quite different from the logarithmic 
laws predicted by the latter. 

Our conclusion is that the SRSB theory seems to be, at the moment, the best picture 
to describe the EA model in finite dimensions greater than the lower critical dimension. 

A still open problem we are planning to study in the future regards the estimate, simu- 
lating much larger lattices, of the order parameter using the static spin glass susceptibility. 
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